Comparing newt and snake phenotype correlation

newt and snake data from Hanifin et al. 2008

snewt <- read.csv("~/Desktop/snake_newt_real_data.csv", header = TRUE)

snewt$mean_res = snewt$mean_toxin+snewt$dif
cor(snewt$mean_toxin, snewt$mean_res)
## [1] 0.7660207
(sites <- data.frame(longitude = snewt$lon, latitude = snewt$lat))
##    longitude latitude
## 1   124.4964  49.9415
## 2   124.5653  49.7453
## 3   125.4494  49.6506
## 4   120.4164  37.2972
## 5   124.2026  41.7558
## 6   121.9018  37.8534
## 7   121.5566  37.0030
## 8   121.7326  40.5401
## 9   120.5733  38.5813
## 10  124.0598  41.2864
## 11  122.3255  37.5630
## 12  121.1893  35.6440
## 13  120.5724  34.7420
## 14  123.3556  39.4096
## 15  123.6314  40.9396
## 16  116.6764  46.7360
## 17  123.3874  44.6282
## 18  124.3873  42.7068
## 19  122.4554  42.1040
## 20  121.8897  44.2164
## 21  123.5642  43.0962
## 22  123.9238  46.1651
## 23  123.9454  46.2729
## 24  124.0168  48.0405
## 25  122.6448  47.1003
## 26  123.8839  46.9770
## 27  122.3490  48.3387
## 28  122.2156  45.6487
just_loc <- data.frame(snewt$lon, snewt$lat)
snewt$lon<- snewt$lon *-1
locations_sf <- st_as_sf(snewt, coords = c("lon", "lat"), crs = 4326)
mapview(locations_sf)
library("rnaturalearth")
library("rnaturalearthdata")
library("tools")
library("maps")
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
world <- ne_countries(scale = "medium", returnclass = "sf")
states <- st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))
states$ID <- toTitleCase(states$ID)

ggplot(data = world) +
  geom_sf() +
  geom_sf(data = states, fill = NA)+
  geom_point(data = snewt, aes(x = lon, y = lat, fill = mean_res), size = 6, shape = 22, ) +
  geom_point(data = snewt, aes(x = lon, y = lat, fill = mean_toxin), size = 4, shape = 21, ) +
  scale_fill_viridis_c()+theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8))+
  coord_sf(xlim = c(-127, -115), ylim = c(34, 50.5), expand = FALSE) 

ggplot(data = world) +
  geom_sf() +
  geom_sf(data = states, fill = NA)+
  #geom_point(data = snewt, aes(x = lon, y = lat, fill = mean_res), size = 6, shape = 22, ) +
  geom_point(data = snewt, aes(x = lon, y = lat, fill = mean_toxin), size = 4, shape = 21, ) +
  scale_fill_viridis_c()+theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8))+
  coord_sf(xlim = c(-127, -115), ylim = c(34, 50.5), expand = FALSE)

ggplot(data = world) +
  geom_sf() +
  geom_sf(data = states, fill = NA)+
  geom_point(data = snewt, aes(x = lon, y = lat, fill = mean_res), size = 6, shape = 22, ) +
  #geom_point(data = snewt, aes(x = lon, y = lat, fill = mean_toxin), size = 4, shape = 21, ) +
  scale_fill_viridis_c()+theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8))+
  coord_sf(xlim = c(-127, -115), ylim = c(34, 50.5), expand = FALSE)

More data from the slim simulation

might be a specific generation or a small set of generations, If I can get the data to be saved and then interpreted by R make a list in a column how to proportion and recreate what I see in the slim map (the different box and their dims)

GA1 experiment values:

GA2 experiment values:

#cor_grid_30 <- list.files(path="~/Desktop/GA1_rep", pattern="_grid_", full.names=TRUE, recursive=FALSE)
#cor_grid_30_G2 <- list.files(path="~/Desktop/GA2_rep", pattern="_grid_", full.names=TRUE, recursive=FALSE)
cor_grid_30 <- list.files(path="~/Desktop/GA1_csv_notperarea", pattern="_grid_", full.names=TRUE, recursive=FALSE)
cor_grid_30_G2 <- list.files(path="~/Desktop/GA2_csv_notperarea", pattern="_grid_", full.names=TRUE, recursive=FALSE)

cor_grid_30_file<-read.csv(cor_grid_30[1], header = TRUE)
replica <- unlist(str_split(unlist(str_split(cor_grid_30[1],"rep_"))[3], "_"))[1]
cor_grid_30_file$rep<-replica
cor_grid_30_file$sim <-"A"
cor_grid_30_file$rep <-unlist(str_split(unlist(str_split(cor_grid_30[1],"rep_"))[3], "_"))[1]
cor_grid_30_file$snake_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30[1],"snake_mu_rate_"))[2], "_"))[1]
cor_grid_30_file$newt_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30[1],"newt_mu_rate_"))[2], "_"))[1]
cor_grid_30_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[1],"snake_mu_effect_sd_"))[2], "_"))[1]
cor_grid_30_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[1],"newt_mu_effect_sd_"))[2], "_"))[1]

cor_grid_30_G2_file<-read.csv(cor_grid_30_G2[1], header = TRUE)
replica <- unlist(str_split(unlist(str_split(cor_grid_30_G2[1],"rep_"))[3], "_"))[1]
cor_grid_30_G2_file$rep<-replica
cor_grid_30_G2_file$sim <-"A"
cor_grid_30_G2_file$rep <-unlist(str_split(unlist(str_split(cor_grid_30_G2[1],"rep_"))[3], "_"))[1]
cor_grid_30_G2_file$snake_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30_G2[1],"snake_mu_rate_"))[2], "_"))[1]
cor_grid_30_G2_file$newt_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30_G2[1],"newt_mu_rate_"))[2], "_"))[1]
cor_grid_30_G2_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30_G2[1],"snake_mu_effect_sd_"))[2], "_"))[1]
cor_grid_30_G2_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30_G2[1],"newt_mu_effect_sd_"))[2], "_"))[1]


oldstand=0
letter=1
for(i in seq(1, length(cor_grid_30_G2), 1)){
  stand <- unlist(str_split(unlist(str_split(cor_grid_30_G2[i],"newt_mu_rate_"))[2], "_"))[1]
  if (stand!=oldstand){
    oldstand=stand
    #letter = letter + 1
  }
  letter = i
  GA2_lit_temp_file <- read.csv(cor_grid_30_G2[i], header = TRUE)
  GA2_lit_temp_file$sim <-LETTERS[letter]
  GA2_lit_temp_file$rep <-unlist(str_split(unlist(str_split(cor_grid_30_G2[i],"rep_"))[3], "_"))[1]
  GA2_lit_temp_file$snake_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30_G2[i],"snake_mu_rate_"))[2], "_"))[1]
GA2_lit_temp_file$newt_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30[i],"newt_mu_rate_"))[2], "_"))[1]
GA2_lit_temp_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[i],"snake_mu_effect_sd_"))[2], "_"))[1]
GA2_lit_temp_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  
  cor_grid_30_G2_file <- rbind(cor_grid_30_G2_file, GA2_lit_temp_file)
 
  
  if (i<=length(cor_grid_30)){
  cor_temp_file <- read.csv(cor_grid_30[i], header = TRUE)
  cor_temp_file$sim <-LETTERS[letter]
  cor_temp_file$rep <-unlist(str_split(unlist(str_split(cor_grid_30[i],"rep_"))[3], "_"))[1]
  cor_temp_file$snake_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30[i],"snake_mu_rate_"))[2], "_"))[1]
cor_temp_file$newt_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30[i],"newt_mu_rate_"))[2], "_"))[1]
cor_temp_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[i],"snake_mu_effect_sd_"))[2], "_"))[1]
cor_temp_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  cor_grid_30_file <- rbind(cor_grid_30_file, cor_temp_file)
  rep <-unlist(str_split(unlist(str_split(cor_grid_30[i],"rep_"))[3], "_"))[1]
  
  snake_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30[i],"snake_mu_rate_"))[2], "_"))[1]
  newt_mu_rate <- unlist(str_split(unlist(str_split(cor_grid_30[i],"newt_mu_rate_"))[2], "_"))[1]
  snake_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[i],"snake_mu_effect_sd_"))[2], "_"))[1]
  newt_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_grid_30[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  print(paste0("Simulation ", LETTERS[letter],": Snake mu-rate & effect sd (", snake_mu_rate,", ", snake_mu_effect_sd, ") Newt mu-rate & effect sd (", newt_mu_rate, ", ", newt_mu_effect_sd, ")",  " rep: ", rep))
  }
  
}
## [1] "Simulation A: Snake mu-rate & effect sd (1.0e-08, 0.005) Newt mu-rate & effect sd (1.0e-08, 0.005) rep: 0"
## [1] "Simulation B: Snake mu-rate & effect sd (1.0e-08, 0.005) Newt mu-rate & effect sd (1.0e-09, 0.05) rep: 0"
## [1] "Simulation C: Snake mu-rate & effect sd (1.0e-08, 0.005) Newt mu-rate & effect sd (1.0e-10, 0.5) rep: 0"
## [1] "Simulation D: Snake mu-rate & effect sd (1.0e-08, 0.005) Newt mu-rate & effect sd (1.0e-11, 5.0) rep: 0"
## [1] "Simulation E: Snake mu-rate & effect sd (1.0e-09, 0.05) Newt mu-rate & effect sd (1.0e-08, 0.005) rep: 0"
## [1] "Simulation F: Snake mu-rate & effect sd (1.0e-09, 0.05) Newt mu-rate & effect sd (1.0e-09, 0.05) rep: 0"
## [1] "Simulation G: Snake mu-rate & effect sd (1.0e-09, 0.05) Newt mu-rate & effect sd (1.0e-10, 0.5) rep: 0"
## [1] "Simulation H: Snake mu-rate & effect sd (1.0e-09, 0.05) Newt mu-rate & effect sd (1.0e-11, 5.0) rep: 0"
## [1] "Simulation I: Snake mu-rate & effect sd (1.0e-10, 0.5) Newt mu-rate & effect sd (1.0e-08, 0.005) rep: 0"
## [1] "Simulation J: Snake mu-rate & effect sd (1.0e-10, 0.5) Newt mu-rate & effect sd (1.0e-09, 0.05) rep: 0"
## [1] "Simulation K: Snake mu-rate & effect sd (1.0e-10, 0.5) Newt mu-rate & effect sd (1.0e-10, 0.5) rep: 0"
## [1] "Simulation L: Snake mu-rate & effect sd (1.0e-10, 0.5) Newt mu-rate & effect sd (1.0e-11, 5.0) rep: 0"
## [1] "Simulation M: Snake mu-rate & effect sd (1.0e-11, 5.0) Newt mu-rate & effect sd (1.0e-08, 0.005) rep: 0"
## [1] "Simulation N: Snake mu-rate & effect sd (1.0e-11, 5.0) Newt mu-rate & effect sd (1.0e-09, 0.05) rep: 0"
## [1] "Simulation O: Snake mu-rate & effect sd (1.0e-11, 5.0) Newt mu-rate & effect sd (1.0e-10, 0.5) rep: 0"
## [1] "Simulation P: Snake mu-rate & effect sd (1.0e-11, 5.0) Newt mu-rate & effect sd (1.0e-11, 5.0) rep: 0"
cor_grid_30_file$snake <- paste0(cor_grid_30_file$snake_mu_rate, "_",cor_grid_30_file$snake_mu_effect_sd)
cor_grid_30_file$newt <- paste0(cor_grid_30_file$newt_mu_rate, "_",cor_grid_30_file$newt_mu_effect_sd)

cor_grid_30_file$newt_snake_rep <- paste0("Newt_",cor_grid_30_file$newt_mu_rate, "_", cor_grid_30_file$newt_mu_effect_sd, "_Snake_",cor_grid_30_file$snake_mu_rate, "_", cor_grid_30_file$snake_mu_effect_sd, "_", cor_grid_30_file$rep)



#GA1_mean = cor_grid_30_file[ , grepl( "mean" , names( cor_grid_30_file ) ) ]


#unique(cor_grid_30_file$newt_snake_rep)

GA1<-cor_grid_30_file[which(cor_grid_30_file$generation>=10000 & cor_grid_30_file$gen<=110000 ),]
#GA2<-cor_grid_30_G2_file[which(cor_grid_30_G2_file$generation>=10000 & cor_grid_30_G2_file$gen<=110000 ),]
#length(GA1$generation)


#GA1<-cor_grid_30_file[which(cor_grid_30_file$generation<=50 & cor_grid_30_file$newt_snake_rep=="Newt_1.0e-08_0.005_Snake_1.0e-09_0.05_2"),]

GA1_mean = GA1[,grepl("generation|mean", names(GA1)) ]

#How to remake the dataframe into gen, newt_value, snake_value, location (x, y)?

GA1_t <- data.frame(t(GA1_mean))

loc_x = 1
loc_y = 5
skip = 1
GA1_t$loc_x = 0
GA1_t$loc_y = 0
GA1_t$species = 0
grid_info <- as.character(rownames(GA1_t))
for (i in 2:nrow(GA1_t)){
   #print(as.numeric(unlist(str_split(i,"_pheno_"))[2]))
  
  GA1_t$species[i] = unlist(str_split(grid_info[i],"_mean"))[1]
  
  if(loc_y==0){
    loc_x = loc_x + 1
    loc_y = 5
  }
  #print(paste0(loc_x, " ", loc_y))
  GA1_t$loc_x[i]= loc_x
GA1_t$loc_y[i] = loc_y
  
  if (skip %% 2 == 0 ){
    loc_y = loc_y - 1
  } 
  skip = skip +1
  
}

GA1_t$names<- rownames(GA1_t)

GA1_mean = cor_grid_30_file[,grepl( "mean" , names(cor_grid_30_file)) ]
GA1_mean$generation <- cor_grid_30_file$generation
GA1_mean$rep <- cor_grid_30_file$rep
GA1_mean$snake_mu_rate <- cor_grid_30_file$snake_mu_rate
GA1_mean$newt_mu_rate <- cor_grid_30_file$newt_mu_rate
GA1_mean$snake_mu_effect_sd <- cor_grid_30_file$snake_mu_effect_sd
GA1_mean$newt_mu_effect_sd <- cor_grid_30_file$newt_mu_effect_sd
GA1_mean$snake <- cor_grid_30_file$snake
GA1_mean$newt <- cor_grid_30_file$newt
GA1_mean$newt_snake_rep <- cor_grid_30_file$newt_snake_rep


data_long_mean <- gather(GA1_mean, condition, measurement, newt_mean_pheno_0:snake_mean_pheno_24, factor_key=TRUE)


data_long_mean$lox <- GA1_t[match(data_long_mean$condition, GA1_t$names),"loc_x"]
data_long_mean$loy <- GA1_t[match(data_long_mean$condition, GA1_t$names),"loc_y"]
data_long_mean$species <- GA1_t[match(data_long_mean$condition, GA1_t$names),"species"]



GA1_max = cor_grid_30_file[,grepl( "max" , names(cor_grid_30_file)) ]
GA1_max$generation <- cor_grid_30_file$generation
GA1_max$rep <- cor_grid_30_file$rep
GA1_max$snake_mu_rate <- cor_grid_30_file$snake_mu_rate
GA1_max$newt_mu_rate <- cor_grid_30_file$newt_mu_rate
GA1_max$snake_mu_effect_sd <- cor_grid_30_file$snake_mu_effect_sd
GA1_max$newt_mu_effect_sd <- cor_grid_30_file$newt_mu_effect_sd
GA1_max$snake <- cor_grid_30_file$snake
GA1_max$newt <- cor_grid_30_file$newt
GA1_max$newt_snake_rep <- cor_grid_30_file$newt_snake_rep

data_long_max <- gather(GA1_max, condition, measurement, newt_max_pheno_0:snake_max_pheno_24, factor_key=TRUE)

data_long_max$lox <- data_long_mean$lox 
data_long_max$loy <- data_long_mean$loy
data_long_max$species <- data_long_mean$species
GA1_onefile<-data_long_mean[which(data_long_mean$newt_snake_rep=="Newt_1.0e-10_0.5_Snake_1.0e-10_0.5_0" & data_long_max$gen==10981 ),]

#GA1_onefile<-data_long_mean[which(data_long_mean$newt_snake_rep==unique(data_long_mean$newt_snake_rep)[16] & data_long_max$gen==29981 ),]

#GA1_onefile<-data_long_max[which(data_long_max$newt_snake_rep=="Newt_1.0e-10_0.5_Snake_1.0e-10_0.5_0" & data_long_max$gen==29051 ),]

ggplot(GA1_onefile)+
 
  geom_point(data = GA1_onefile[which(GA1_onefile$species=="snake"),], aes(x = lox, y = loy, fill = measurement), size = 10, shape = 22, ) +
  geom_point(data = GA1_onefile[which(GA1_onefile$species=="newt"),], aes(x = lox, y = loy, fill = measurement), size = 7, shape = 21, ) +
  scale_fill_viridis_c()+
  theme_bw()

unique(data_long_mean$newt_snake_rep)
##  [1] "Newt_1.0e-08_0.005_Snake_1.0e-08_0.005_0"
##  [2] "Newt_1.0e-09_0.05_Snake_1.0e-08_0.005_0" 
##  [3] "Newt_1.0e-10_0.5_Snake_1.0e-08_0.005_0"  
##  [4] "Newt_1.0e-11_5.0_Snake_1.0e-08_0.005_0"  
##  [5] "Newt_1.0e-08_0.005_Snake_1.0e-09_0.05_0" 
##  [6] "Newt_1.0e-09_0.05_Snake_1.0e-09_0.05_0"  
##  [7] "Newt_1.0e-10_0.5_Snake_1.0e-09_0.05_0"   
##  [8] "Newt_1.0e-11_5.0_Snake_1.0e-09_0.05_0"   
##  [9] "Newt_1.0e-08_0.005_Snake_1.0e-10_0.5_0"  
## [10] "Newt_1.0e-09_0.05_Snake_1.0e-10_0.5_0"   
## [11] "Newt_1.0e-10_0.5_Snake_1.0e-10_0.5_0"    
## [12] "Newt_1.0e-11_5.0_Snake_1.0e-10_0.5_0"    
## [13] "Newt_1.0e-08_0.005_Snake_1.0e-11_5.0_0"  
## [14] "Newt_1.0e-09_0.05_Snake_1.0e-11_5.0_0"   
## [15] "Newt_1.0e-10_0.5_Snake_1.0e-11_5.0_0"    
## [16] "Newt_1.0e-11_5.0_Snake_1.0e-11_5.0_0"
GA1_onefile_2<-data_long_mean[which(data_long_mean$newt_snake_rep=="Newt_1.0e-08_0.005_Snake_1.0e-08_0.005_0" & data_long_max$gen==10981 ),]



ggplot(GA1_onefile_2)+
 
  geom_point(data = GA1_onefile_2[which(GA1_onefile_2$species=="snake"),], aes(x = lox, y = loy, fill = measurement), size = 10, shape = 22, ) +
  geom_point(data = GA1_onefile_2[which(GA1_onefile_2$species=="newt"),], aes(x = lox, y = loy, fill = measurement), size = 7, shape = 21, ) +
  scale_fill_viridis_c()+
  theme_bw()

snake_info<-data_long_mean[which(data_long_mean$species=="snake"),]
newt_info<-data_long_mean[which(data_long_mean$species=="newt"),]

names(newt_info)[names(newt_info) == 'measurement'] <- 'newt_toxin'
newt_info$snake_res<-snake_info$measurement
newt_info$dif<-newt_info$snake_res-newt_info$newt_toxin
pattern<-"Newt_1.0e-08_0.005_Snake_1.0e-08_0.005_0"
scatter_pheno<-newt_info[which(newt_info$newt_snake_rep==pattern & newt_info$generation>14000 & newt_info$generation<15000 ),]
#scatter_pheno<-newt_info[which(newt_info$newt_snake_rep==pattern & newt_info$generation==20981 ),]
scatter_pheno[scatter_pheno == 0] <- NA
scatter_pheno[scatter_pheno == 1] <- NA
ggplot(scatter_pheno)+
 ggtitle(paste0(" pattern ", pattern, " From ", min(scatter_pheno$generation), " gen to ", max(scatter_pheno$generation), " gen"))+
  geom_point(aes(y = newt_toxin, x = snake_res, color=generation), size = 3) 
## Warning: Removed 350 rows containing missing values (`geom_point()`).

  theme_bw()
## List of 97
##  $ line                      :List of 6
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ lineend      : chr "butt"
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ rect                      :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ text                      :List of 11
##   ..$ family       : chr ""
##   ..$ face         : chr "plain"
##   ..$ colour       : chr "black"
##   ..$ size         : num 11
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 0
##   ..$ lineheight   : num 0.9
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ title                     : NULL
##  $ aspect.ratio              : NULL
##  $ axis.title                : NULL
##  $ axis.title.x              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.75points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.top          :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.75points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.bottom       : NULL
##  $ axis.title.y              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.75points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.y.left         : NULL
##  $ axis.title.y.right        :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.75points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text                 :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey30"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.2points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.top           :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.2points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.bottom        : NULL
##  $ axis.text.y               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 1
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.2points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.y.left          : NULL
##  $ axis.text.y.right         :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.2points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.ticks                :List of 6
##   ..$ colour       : chr "grey20"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ axis.ticks.x              : NULL
##  $ axis.ticks.x.top          : NULL
##  $ axis.ticks.x.bottom       : NULL
##  $ axis.ticks.y              : NULL
##  $ axis.ticks.y.left         : NULL
##  $ axis.ticks.y.right        : NULL
##  $ axis.ticks.length         : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  $ axis.ticks.length.x       : NULL
##  $ axis.ticks.length.x.top   : NULL
##  $ axis.ticks.length.x.bottom: NULL
##  $ axis.ticks.length.y       : NULL
##  $ axis.ticks.length.y.left  : NULL
##  $ axis.ticks.length.y.right : NULL
##  $ axis.line                 : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.line.x               : NULL
##  $ axis.line.x.top           : NULL
##  $ axis.line.x.bottom        : NULL
##  $ axis.line.y               : NULL
##  $ axis.line.y.left          : NULL
##  $ axis.line.y.right         : NULL
##  $ legend.background         :List of 5
##   ..$ fill         : NULL
##   ..$ colour       : logi NA
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ legend.margin             : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing            : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing.x          : NULL
##  $ legend.spacing.y          : NULL
##  $ legend.key                :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : logi NA
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ legend.key.size           : 'simpleUnit' num 1.2lines
##   ..- attr(*, "unit")= int 3
##  $ legend.key.height         : NULL
##  $ legend.key.width          : NULL
##  $ legend.text               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.text.align         : NULL
##  $ legend.title              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.title.align        : NULL
##  $ legend.position           : chr "right"
##  $ legend.direction          : NULL
##  $ legend.justification      : chr "center"
##  $ legend.box                : NULL
##  $ legend.box.just           : NULL
##  $ legend.box.margin         : 'margin' num [1:4] 0cm 0cm 0cm 0cm
##   ..- attr(*, "unit")= int 1
##  $ legend.box.background     : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.box.spacing        : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##  $ panel.background          :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : logi NA
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ panel.border              :List of 5
##   ..$ fill         : logi NA
##   ..$ colour       : chr "grey20"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ panel.spacing             : 'simpleUnit' num 5.5points
##   ..- attr(*, "unit")= int 8
##  $ panel.spacing.x           : NULL
##  $ panel.spacing.y           : NULL
##  $ panel.grid                :List of 6
##   ..$ colour       : chr "grey92"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ panel.grid.major          : NULL
##  $ panel.grid.minor          :List of 6
##   ..$ colour       : NULL
##   ..$ linewidth    : 'rel' num 0.5
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ panel.grid.major.x        : NULL
##  $ panel.grid.major.y        : NULL
##  $ panel.grid.minor.x        : NULL
##  $ panel.grid.minor.y        : NULL
##  $ panel.ontop               : logi FALSE
##  $ plot.background           :List of 5
##   ..$ fill         : NULL
##   ..$ colour       : chr "white"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ plot.title                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 1.2
##   ..$ hjust        : num 0
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 5.5points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.title.position       : chr "panel"
##  $ plot.subtitle             :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 5.5points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.caption              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : num 1
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 5.5points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.caption.position     : chr "panel"
##  $ plot.tag                  :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 1.2
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.tag.position         : chr "topleft"
##  $ plot.margin               : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
##   ..- attr(*, "unit")= int 8
##  $ strip.background          :List of 5
##   ..$ fill         : chr "grey85"
##   ..$ colour       : chr "grey20"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ strip.background.x        : NULL
##  $ strip.background.y        : NULL
##  $ strip.clip                : chr "inherit"
##  $ strip.placement           : chr "inside"
##  $ strip.text                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey10"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 4.4points 4.4points 4.4points 4.4points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.text.x              : NULL
##  $ strip.text.x.bottom       : NULL
##  $ strip.text.x.top          : NULL
##  $ strip.text.y              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.text.y.left         :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.text.y.right        : NULL
##  $ strip.switch.pad.grid     : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  $ strip.switch.pad.wrap     : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi TRUE
##  - attr(*, "validate")= logi TRUE
cor_files <- list.files(path="~/Desktop/csv_GA1_long", pattern="_cor_", full.names=TRUE, recursive=FALSE)
cor_GA2_files <- list.files(path="~/Desktop/GA2_rep", pattern="_cor_", full.names=TRUE, recursive=FALSE)
lit_30 <- list.files(path="~/Desktop/csv_GA1_long", pattern="_lit_", full.names=TRUE, recursive=FALSE)
lit_30_G2 <- list.files(path="~/Desktop/GA2_rep", pattern="_lit_", full.names=TRUE, recursive=FALSE)

lit_30_file<- data.frame()
lit_30_G2_file <- data.frame()


number=1
#cor_file <- read.csv(cor_files[number], header = TRUE)
#cor_file$sim <-"A"
#cor_GA2_file <- read.csv(cor_GA2_files[number], header = TRUE)
#cor_GA2_file$sim <-"A"

cor_file<- data.frame()
cor_GA2_file<- data.frame()
oldstand=0
letter=0
newt_unique=0
newt_numb = 0
for(i in seq(1, length(cor_GA2_files), 1)){
  stand <- unlist(str_split(unlist(str_split(cor_GA2_files[i],"newt_mu_rate_"))[2], "_"))[1]
  up <- unlist(str_split(unlist(str_split(cor_GA2_files[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  phrase = paste0(stand, "_", up)
  if (stand!=oldstand){
    oldstand=stand
    letter = letter + 1
  }

    if (phrase=="1.0e-08_0.05" ){
    newt_numb =  1
   }
    if (phrase=="1.0e-09_0.158114"){
    newt_numb =  2
   }
      if (phrase=="1.0e-10_0.5"){
    newt_numb =  3
   }
    if (phrase=="1.0e-11_1.58114"){
    newt_numb =  4
   }
  
      if (phrase=="1.0e-12_5.0"){
    newt_numb =  5
   }

  

  GA2_cor_temp_file <- read.csv(cor_GA2_files[i], header = TRUE)
  GA2_cor_temp_file$sim <-LETTERS[letter]
  GA2_cor_temp_file$size <-newt_numb
  GA2_cor_temp_file$rep <-unlist(str_split(unlist(str_split(cor_GA2_files[i],"rep_"))[3], "_"))[1]
  GA2_cor_temp_file$snake_mu_rate <- unlist(str_split(unlist(str_split(cor_GA2_files[i],"snake_mu_rate_"))[2], "_"))[1]
  GA2_cor_temp_file$newt_mu_rate <- unlist(str_split(unlist(str_split(cor_GA2_files[i],"newt_mu_rate_"))[2], "_"))[1]
  GA2_cor_temp_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_GA2_files[i],"snake_mu_effect_sd_"))[2], "_"))[1]
  GA2_cor_temp_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_GA2_files[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  
  cor_GA2_file <- rbind(cor_GA2_file, GA2_cor_temp_file)
  
  GA2_lit_temp_file <- read.csv(lit_30_G2[i], header = TRUE)
  GA2_lit_temp_file$sim <-LETTERS[letter]
  GA2_lit_temp_file$size <-newt_numb
  GA2_lit_temp_file$rep <-unlist(str_split(unlist(str_split(lit_30_G2[i],"rep_"))[3], "_"))[1]
  GA2_lit_temp_file$snake_mu_rate <- unlist(str_split(unlist(str_split(lit_30_G2[i],"snake_mu_rate_"))[2], "_"))[1]
  GA2_lit_temp_file$newt_mu_rate <- unlist(str_split(unlist(str_split(lit_30_G2[i],"newt_mu_rate_"))[2], "_"))[1]
  GA2_lit_temp_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(lit_30_G2[i],"snake_mu_effect_sd_"))[2], "_"))[1]
  GA2_lit_temp_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(lit_30_G2[i],"newt_mu_effect_sd_"))[2], "_"))[1]

  
  lit_30_G2_file <- rbind(lit_30_G2_file, GA2_lit_temp_file)
  

  if (i<=length(cor_files)){
  stand <- unlist(str_split(unlist(str_split(cor_files[i],"newt_mu_rate_"))[2], "_"))[1]
  up <- unlist(str_split(unlist(str_split(cor_files[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  phrase = paste0(stand, "_", up)

   if (phrase=="1.0e-08_0.005"){
    newt_numb =  1
   }
    if (phrase=="1.0e-09_0.05"){
    newt_numb =  2
   }
      if (phrase=="1.0e-10_0.5"){
    newt_numb =  3
   }
    if (phrase=="1.0e-11_5.0"){
    newt_numb =  4
   }
  cor_temp_file <- read.csv(cor_files[i], header = TRUE)
  cor_temp_file$sim <-LETTERS[letter]
  cor_temp_file$size <-newt_numb
  cor_temp_file$rep <-unlist(str_split(unlist(str_split(cor_files[i],"rep_"))[3], "_"))[1]
  cor_temp_file$snake_mu_rate <- unlist(str_split(unlist(str_split(cor_files[i],"snake_mu_rate_"))[2], "_"))[1]
  cor_temp_file$newt_mu_rate <- unlist(str_split(unlist(str_split(cor_files[i],"newt_mu_rate_"))[2], "_"))[1]
  cor_temp_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_files[i],"snake_mu_effect_sd_"))[2], "_"))[1]
  cor_temp_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(cor_files[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  cor_file <- rbind(cor_file, cor_temp_file)
  
  GA1_temp_file <- read.csv(lit_30[i], header = TRUE)
  GA1_temp_file$sim <-LETTERS[letter]
  GA1_temp_file$size <-newt_numb
  GA1_temp_file$rep <-unlist(str_split(unlist(str_split(lit_30[i],"rep_"))[3], "_"))[1]
  GA1_temp_file$snake_mu_rate <- unlist(str_split(unlist(str_split(lit_30[i],"snake_mu_rate_"))[2], "_"))[1]
  GA1_temp_file$newt_mu_rate <- unlist(str_split(unlist(str_split(lit_30[i],"newt_mu_rate_"))[2], "_"))[1]
  GA1_temp_file$snake_mu_effect_sd <- unlist(str_split(unlist(str_split(lit_30[i],"snake_mu_effect_sd_"))[2], "_"))[1]
  GA1_temp_file$newt_mu_effect_sd <- unlist(str_split(unlist(str_split(lit_30[i],"newt_mu_effect_sd_"))[2], "_"))[1]
  lit_30_file <- rbind(lit_30_file, GA1_temp_file)
  
  #print(paste0("Simulation ", LETTERS[letter],": Snake mu-rate & effect sd (", cor_temp_file$snake_mu_rate[i],", ", cor_temp_file$snake_mu_effect_sd[i], ") Newt mu-rate & effect sd (", cor_temp_file$newt_mu_rate[i], ", ",cor_temp_file$ newt_mu_effect_sd[i], ")" ))
  }
  
}
lit_30_file$dif <- lit_30_file$Snake_mean_Pheno-lit_30_file$Newt_mean_Pheno
lit_30_G2_file$dif <- lit_30_G2_file$Snake_mean_Pheno-lit_30_G2_file$Newt_mean_Pheno
boxplot_cors_by_gen <- function(file, title, subt, x_label, y_label, min_gen, max_gen){
  real_snewt_cor = cor(snewt$mean_toxin, snewt$mean_res)
  file$newt_snake_rep <- paste0(file$newt_mu_rate, "_", file$newt_mu_effect_sd, "_", file$snake_mu_rate, "_", file$snake_mu_effect_sd, "_", file$rep)
  file$newt_snake <- paste0(file$newt_mu_rate, "_", file$newt_mu_effect_sd, "_", file$snake_mu_rate, "_", file$snake_mu_effect_sd)
  unique(file$sim)
  plot <- ggplot(file[which(file$gen>=min_gen & file$gen<=max_gen),]) +
  ylim(-1.1, 1.1)+
  ggtitle(title, subtitle = paste0(subt, " From ", min_gen, " gen to ", max_gen, " gen"))+
  xlab(x_label) + ylab(y_label)+
  geom_boxplot(aes(x=factor(rep), y=as.numeric(mean_newt_pheno_By_mean_snake_pheno), color = "CA"))+
  geom_hline(yintercept=real_snewt_cor, linetype="dashed")+
    geom_hline(yintercept=0, linetype="solid")+
  #geom_boxplot(aes(x=factor("CB"), y=as.numeric(num_newts_By_num_snakes), color = "CB"))+
  #geom_boxplot(aes(x=factor("CC"), y=as.numeric(sum_newt_pheno_By_num_snake), color = "CC"))+
  #geom_boxplot(aes(x=factor("CD"), y=as.numeric(sum_snake_pheno_By_num_newt), color = "CD"))+
  theme(plot.title = element_text(size=5)) +
  facet_wrap(~newt_snake)+
  theme_bw() +  theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8), legend.position="none",
                      strip.background = element_blank(), strip.text.x = element_blank(),
                      title =element_text(size=8, face='bold'), axis.title.x = element_blank(),
                      axis.title.y = element_blank())


  return(plot)
}

var_cor <- function(file, cfile, title, subt, x_label, y_label, min_gen, max_gen){
  real_snewt_cor = cor(snewt$mean_toxin, snewt$mean_res)
  real_snewt_var = var(snewt$dif)
  cfile$newt_snake_rep <- paste0("Newt_", cfile$newt_mu_rate, "_", cfile$newt_mu_effect_sd, "_Snake_", cfile$snake_mu_rate, "_", cfile$snake_mu_effect_sd, "_", cfile$rep)
  var_cor_table=data.frame()

  
  for(i in 1:length(unique(file$newt_snake_rep))){
    pattern<-unique(file$newt_snake_rep)[i]
    temp_file<-file[which(file$gen>=min_gen & file$gen<=max_gen & file$newt_snake_rep==pattern),]
    temp_cor_file<-cfile[which(cfile$gen>=min_gen & cfile$gen<=max_gen & cfile$newt_snake_rep==pattern),]
    
    var_cor_table_temp<- data.frame("pattern"=pattern, "var"= var(temp_file$dif),"cor" = cor(temp_file$newt_toxin, temp_file$snake_res), "mean_corr" = mean(temp_cor_file$mean_newt_pheno_By_mean_snake_pheno))
    var_cor_table <- rbind(var_cor_table, var_cor_table_temp)
  }
 
  var_cor_table_temp<- data.frame("pattern"="real", "var"= real_snewt_var, "cor" = real_snewt_cor, "mean_corr" = real_snewt_cor)
  var_cor_table <- rbind(var_cor_table, var_cor_table_temp)
  
  plot <- ggplot(var_cor_table) +
  ggtitle(title, subtitle = paste0(subt, " From ", min_gen, " gen to ", max_gen, " gen"))+
  xlab(x_label) + ylab(y_label)+
  geom_point(aes(x=as.numeric(var), y=as.numeric(cor), color = pattern), shape=17, size=3)+
  geom_point(aes(x=as.numeric(var), y=as.numeric(mean_corr), color = pattern), size=3)+
  theme(plot.title = element_text(size=5)) +
  theme_bw() +  theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8))


  return(plot)
}



trends_by_gen <- function(file, cfile, title, subt, x_label, x_label_cor, y_label, min_gen, max_gen, pattern){
 file$newt_snake_rep <- paste0(file$newt_mu_rate, "_", file$newt_mu_effect_sd, "_", file$snake_mu_rate, "_", file$snake_mu_effect_sd, "_", file$rep)
  cfile$newt_snake_rep <- paste0(cfile$newt_mu_rate, "_", cfile$newt_mu_effect_sd, "_", cfile$snake_mu_rate, "_", cfile$snake_mu_effect_sd, "_", cfile$rep)
  plot <- ggplot(file[which(file$gen>=min_gen & file$gen<=max_gen & file$newt_snake_rep==pattern),]) +
    #geom_line(aes(generation, dif, group=rep, col="dif"))+
      ggtitle(title, subtitle = paste0(subt," pattern ", pattern, " From ", min_gen, " gen to ", max_gen, " gen"))+
    geom_line(aes(generation, Snake_mean_Pheno, group=rep, color="Snake_mean_Pheno"),)+
    geom_line(aes(generation, Newt_mean_Pheno, group=rep, color="Newt_mean_Pheno"))+
    #ylim(0,10)+
    #xlab(x_label) + ylab(y_label)+
    scale_color_manual(name = "Type",
values = c( "Snake_mean_Pheno" = "steelblue", "Newt_mean_Pheno" = "darkred"),
labels = c("Snake", "Newt"))+
    theme_bw() +theme(legend.position="none",
                      strip.background = element_blank(), strip.text.x = element_blank(),
                      title =element_text(size=8, face='bold'), axis.title.x = element_blank(),
                      axis.title.y = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank(),
                      plot.margin=unit(c(1,1,+0.3,1), "cm"))
  
  cor_plot <- ggplot(cfile[which(cfile$gen>=min_gen & cfile$gen<=max_gen & cfile$newt_snake_rep==pattern),]) +
    geom_line(aes(generation, mean_newt_pheno_By_mean_snake_pheno, col="CA"))+
    #ggtitle(title, subtitle = subt) +
    ylim(-1.1, 1.1)+
    #xlab(x_label_cor) + ylab(y_label)+
    theme_bw() +theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8), legend.position="none",
                      strip.background = element_blank(), strip.text.x = element_blank(),
                      title =element_text(size=8, face='bold'), axis.title.x = element_blank(),
                      axis.title.y = element_blank(),
                      plot.margin=unit(c(-0.5,1,1,1), "cm"))
  
  file<-file[which(file$gen>=min_gen & file$gen<=max_gen & file$newt_snake_rep==pattern),]
  cfile<-cfile[which(cfile$gen>=min_gen & cfile$gen<=max_gen & cfile$newt_snake_rep==pattern),]
  toDelete <- seq(0, nrow(cfile), 2)
  cfile_short <-  cfile[-toDelete, ]
  figure <- ggarrange(plot, cor_plot + font("x.text", size = 10),heights = c(2, 1.5),ncol = 1, nrow = 2, align = "v")
    print(paste0("pattern ", pattern))
  print(paste0("Cor between average snake pheno and local cor ", cor(file$Snake_mean_Pheno, cfile_short$mean_newt_pheno_By_mean_snake_pheno)))
  print(paste0("Cor between average newt pheno and local cor ", cor(file$Newt_mean_Pheno, cfile_short$mean_newt_pheno_By_mean_snake_pheno)))
  print(paste0("Cor between average dif pheno and local cor ", cor(file$dif, cfile_short$mean_newt_pheno_By_mean_snake_pheno)))
  return(figure)
}


four_types_boxplot_cors_by_gen <- function(file,title, subt, x_label, y_label, min_gen, max_gen){
  file$newt_snake <- paste0(file$newt_mu_rate, "_", file$newt_mu_effect_sd, "_", file$snake_mu_rate, "_", file$snake_mu_effect_sd)
  plot <- ggplot(file[which(file$gen>=min_gen & file$gen<=max_gen),]) +
  ylim(-1.1, 1.1)+
  ggtitle(title, subtitle = paste0(subt, " From ", min_gen, " gen to ", max_gen, " gen"))+
  xlab(x_label) + ylab(y_label)+
  geom_boxplot(aes(x=factor("CA"), y=as.numeric(mean_newt_pheno_By_mean_snake_pheno), color = "CA"))+
  geom_boxplot(aes(x=factor("CB"), y=as.numeric(num_newts_By_num_snakes), color = "CB"))+
  geom_boxplot(aes(x=factor("CC"), y=as.numeric(sum_newt_pheno_By_num_snake), color = "CC"))+
  geom_boxplot(aes(x=factor("CD"), y=as.numeric(sum_snake_pheno_By_num_newt), color = "CD"))+
  geom_hline(yintercept=0, linetype="solid")+
  theme(plot.title = element_text(size=5)) +
  facet_wrap(~newt_snake)+
  theme_bw() +  theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8), legend.position="none",
                      strip.background = element_blank(), strip.text.x = element_blank(),
                      title =element_text(size=8, face='bold'), axis.title.x = element_blank(),
                      axis.title.y = element_blank())


  return(plot)
}
boxplot_cors_by_gen(cor_file, "Simulation Correlation Comparied to Real Correlation", "GA1", "newt/snake GAs&rep", "Cor", 4000, 5000)

four_types_boxplot_cors_by_gen(cor_file, "Simulation four Correlation", "GA1", "newt/snake GAs&rep", "Cor", 4000, 5000)

boxplot_cors_by_gen(cor_GA2_file, "Simulation Correlation Comparied to Real Correlation", "GA2", "newt/snake GAs&rep", "Cor", 4000, 5000)

trends_by_gen(lit_30_file, cor_file, "test", "GA1", "newt/snake GAs&rep","other", "Cor", min_gen=0, max_gen=5000, pattern="1.0e-10_0.5_1.0e-10_0.5_0")
## [1] "pattern 1.0e-10_0.5_1.0e-10_0.5_0"
## [1] "Cor between average snake pheno and local cor 0.257236703085727"
## [1] "Cor between average newt pheno and local cor 0.0460065772310917"
## [1] "Cor between average dif pheno and local cor 0.0735100676282853"

trends_by_gen(lit_30_file, cor_file, "test", "GA1", "newt/snake GAs&rep","other", "Cor", min_gen=10000, max_gen=15000, pattern="1.0e-10_0.5_1.0e-10_0.5_0")
## [1] "pattern 1.0e-10_0.5_1.0e-10_0.5_0"
## [1] "Cor between average snake pheno and local cor -0.0244914135160876"
## [1] "Cor between average newt pheno and local cor 0.331733410464774"
## [1] "Cor between average dif pheno and local cor -0.208485161185655"

trends_by_gen(lit_30_file, cor_file, "test", "GA1", "newt/snake GAs&rep","other", "Cor", min_gen=0, max_gen=5000, pattern="1.0e-08_0.005_1.0e-08_0.005_0")
## [1] "pattern 1.0e-08_0.005_1.0e-08_0.005_0"
## [1] "Cor between average snake pheno and local cor 0.103828794302331"
## [1] "Cor between average newt pheno and local cor 0.299141985390285"
## [1] "Cor between average dif pheno and local cor -0.372704622259776"

trends_by_gen(lit_30_file, cor_file, "test", "GA1", "newt/snake GAs&rep","other", "Cor", min_gen=0, max_gen=15000, pattern="1.0e-09_0.05_1.0e-08_0.005_0")
## [1] "pattern 1.0e-09_0.05_1.0e-08_0.005_0"
## [1] "Cor between average snake pheno and local cor 0.0153139143652654"
## [1] "Cor between average newt pheno and local cor -0.0284561836126627"
## [1] "Cor between average dif pheno and local cor 0.0280010323520237"

trends_by_gen(lit_30_file, cor_file, "test", "GA1", "newt/snake GAs&rep","other", "Cor", min_gen=0, max_gen=15000, pattern="1.0e-08_0.005_1.0e-09_0.05_0")
## [1] "pattern 1.0e-08_0.005_1.0e-09_0.05_0"
## [1] "Cor between average snake pheno and local cor 0.0298135650824854"
## [1] "Cor between average newt pheno and local cor 0.0977110523590969"
## [1] "Cor between average dif pheno and local cor -0.00538798771324545"

trends_by_gen(lit_30_G2_file, cor_GA2_file, "test", "GA1", "newt/snake GAs&rep","other", "Cor", min_gen=0, max_gen=15000, pattern="1.0e-08_0.05_1.0e-08_0.05_0")
## [1] "pattern 1.0e-08_0.05_1.0e-08_0.05_0"
## [1] "Cor between average snake pheno and local cor 0.117586087993649"
## [1] "Cor between average newt pheno and local cor 0.0720273011182107"
## [1] "Cor between average dif pheno and local cor 0.0342810647748496"

snake_info<-data_long_mean[which(data_long_mean$species=="snake"),]
newt_info<-data_long_mean[which(data_long_mean$species=="newt"),]

names(newt_info)[names(newt_info) == 'measurement'] <- 'newt_toxin'
newt_info$snake_res<-snake_info$measurement
newt_info$dif<-newt_info$snake_res-newt_info$newt_toxin


var_cor(newt_info, cor_file, "var by cor", "mini", "var", "cor", 10000, 15000)